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АВ5ТКАСТ 


The development of microelectronics has brought into being 
large capacity digital memories in a small package. In the 
foreseeable future even more advances can be seen in this 
trend. Therefore the use of digital computers in control sys- 
tems will play an even larger role than today. 

This work involves a fourth order system to simulate the 
control and dynamics of a missile. Proportional navigation 
is used as the guidance method. Studied are the effects of 
applying different controls which are considered best from a 
computer study and the effects of applying digital filtering 
methods. Although these studies were applied to a specific 
problem, an attempt is made to keep the discussion general 


in order that the methods may be considered for other problems. 
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СНАРТЕК 1 
INTRODUCTION 

If it is desired to hit a target with a missile, some form of 
guidance must be employed, This guidance would be selected 
weighing the desired objectives and an estimate of the target 
capabilities against the amount of money and missile space 
available, As suggested in the introduction, the advances in 
computer size and capabilities is, and will in the future, con- 
tinue to change the weighting of these factors. Thus in the 
foreseeable future it may be possible to design a model of the 
system which will be in sufficient detail to predict the respon- 
se of the real system to a high degree of accuracy. Using the 
micro-computers, this model could then become a part of the 
control system in the missile. Such a control system might be 
described in the block diagram in Figure 1-1. 

In Figure 1-1 the missile measures and reacts to a signal 
which depends on the type of guidance employed. Sucha 
signal might depend on the rate of change in the missile's 
radar seeker head which is tracking the line of sight from the 
missile to target. The missile will react to this signal, pro- 
ducing through its dynamics a change in the direction of its 
velocity vector, or the missile will take some other action 
called for by the particular guidance law being used. The 
model receives the measured signals and considers the 
measurable states from the missile dynamics. Since all the 
states of the model are measurable and the digital model in 


itself generates relatively little noise, the model becomes a 
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Figure lel. Block Diagram of a Dynamic System. 
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state predictor or estimator. 

Since all the states of the model are available, a control 
law considering all the states can now be designed. This 
control law is applied in the controller to guide the missile, 

If the desired objectives of the system can justify the 
above approach, the system performance of the modified system 
should be an improvement over most conventional guidance 


Systems. 
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CHAPTER II 
PROPORTIONAL NAVIGATION 
Adler in ЕЛ defined proportional navigation as, 


A course in which the rate of change of missile heading is 
directly proportional to the rate of rotation of the line-of- 
sight from the missile to the target. 


y= Kg (1) 
where 

y is the angular rate of change of the missile velocity vector. 

gis the angular rate of change of the line-of-sight. 

K is a constant, typically between three and five. 

To gain a better understanding of proportional navigation, 
some of the simpler forms of line-of-sight guidance will be 
examined. Pursuit course (sometimes called pure pursuit) always 
aims the missile directly at the target along the line-of-sight. 
Ihis method will always end in a tail chase, even if the missile 
is launched head-on with the target. High missile acceler- 
ations are required. It can be shown from the equations of 
motion that if the missile to target speed ratiop exceeds two, 
the final missile turning rate will approach infinity. Thus 
pursuit course guidance may not be very satisfactory although 
itis very simple to implement. 

The next step is to lead the line-of-sight by an angle 
which is a function of the target velocity. This type is called 
constant bearing course, and is achieved by aligning the 
relative missile to target velocity vector with the line-of- 
sight. That is to say the line-of-sight maintains a constant 


direction in space. It can be shown that this method can 
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cause a missile to follow a target to a collision even if the 
target is maneuvering. This method however requires an in- 
stantaneous correction to each change in the line-of-sight, 

Consider the constant bearing trajectory shown in 
Figure 2-1. In the figure the line-of-sight is shown for 
several samples made by the radar. As long as the missile 
remains on the collision course, the angular rate of change of 
the line-of-sight will remain zero and no control action is 
required. Any rotation of the line-of-sight, indicating a de- 
parture from the collision course will be detected by the 
missile's radar. The missile using proportional navigation 
will turn at a rate proportional to this rotation in a direction 
to reduce the line-of-sight rate and return to a constant 
bearing collision course. 

In the Figure 2-2 above, the lead angle 8is defined: 

В= о-у 2) 
This angle (8) is independent of any reference used to define 
с andy. If, as in the discussion above, we could achieve a 
constant bearing course, Bwould remain a constant, and 
therefore Bwould remain Zero, 

Thus for a constant bearing course, 

B=a-y=0 (3) 
Since the missile cannot react instantaneously, Bis not 
always zero. Thus for any instant of time (1) may be written 
as, 

Ко- В+ у (4) 
If Bis equal to zero, the missile is on a constant bearing 


collision course. Therefore the object of the control system 
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16 


 >-” «ма» Dee — < w cara ο «ναι, а» u a 


тч». э ermg dë سم‎ == ~ а 


LINE-OF-SIGHT SEEKER 
HEAD 


RADAR | o с Ко У 
DELAY K | DYNAMICS 


νὰ 


Figure 2-3. Missile Guidance Inter Loop. 
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Figure 2-4. Missile Guidance Outer Loop. 
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in proportional navigation is to minimize Band if g becomes a 
constant or Zero to drive B to Zero. 
Guidance and Control 

A typical homing missile inter guidance loop is shown in 
Figure 2-3. The following discussion will be limited to only 
coplanar motion of the missile and target. A constant missile 
and target speed is also assumed throughout the rest of this 
paper. 

In Figure 2-3, the radar tracks the target using a servo 
loop to keep the antenna on the line-of-sight. A delayed sig- 
nal which is a measure ої с, is obtained from the radar system. 
A voltage which is proportional to the rate of change of this 
Signal is obtained. This voltage is compared with the output 
of the dynamic system y, and the error signal is obtained. 
This error signal is sent either to the hydraulic system to 
rotate the control surfaces or to a hydraulic/valve system 
which controls side thrust jets, depending on the type of 
missile being used. 

Kinematics 

As might be suggested by the name inter loop, there is 
also an outer loop as shown in Figure 2~4. In the Figure, 
the guidance and control block is shown as in Figure 2-3. 

The effect of y on the line-of-sight and the effect of the tar- 
get motion on the line-of-sight are represented in the missile- 
target kinematics block. 

This outer loop is not as accessible to measurement 
and control, but the understanding of its effects is necessary 


in order to perform a simulation of the entire system. From 
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the geometry of Figure 2-2, the following equations may be 


written: 
R - V, COS(o) - X, COS(B) (5) 
Ro = -Vt SIN(o) + Ма SIN(B) (6) 
where 


ү; = Target speed (assumed Constant) 

Vin = Missile speed (assumed Constant) 

B=o- y 

R = the line-of-sight distance 

R = the rate of change of the line-of-sight. 
Then equations (5) and (6) may be rewritten: 

R = V, COS(o) - Vm COS(o - y) (7) 

Rg 7 -V4 SIN(o) + Vm SIN(o - y) (8) 
Differentiating (8) with respect to time: 

Ro + Rç= -\, СО5(а) в + Ма СО8(а- У) (в- 7) 


= -Ro - Vy, COS(o - >) у (9) 
Therefore 
2RO + RO = -V,, COS(o - >) Z (10) 
let 
R = -R/T 
and 
V, =R} 


where T is the time to go to the target. Then equation (10) 
may be represented by: 
(-2R/T + RS) o = -V COS(0 - y) y (11) 
thus 
o/y = - V, COS(0 - y) (12) 
R/T (TS - 2) 
Where УК СО5(0- y) is called the effective navigation constant. 


R/T 
19 


Li! 


Since V4 is constant and variations in COS(o - y) and V, are 
small, this term is assumed constant and values of four to 


Six are considered best. 
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CHAPTER III 
FILTER AND CONTROL 
Plant Description 
A linear, time-invariant dynamic system is described in 
the flow graph, Figure 3-1. 
The transfer function of Figure 3-1 may be stated: 
Xout (5) Ы Re) ENEE 


pum NEC и пл a c DE 
Xin (8) 51 +а 5 +... + 45 та) (1) 


Equation (1) may be written in one differentia1 equation using 
matrix notation. 

X=FX+DU, C=BX (2) 
Where 

X is the state vector (n x 1) 

Ú is a scalar of the system inputs and controls 

F is a matrix of constants. For the system described in 
Figure 3-1, F is defined in equation (3). 

0 1 0 š š 0 

0 0 1 : : 0 
БӘС < & (3) 

0 (ena . 1 0 

0 0 0 : 0 1 


= = о о о =а 
τ 22 n 
D is a matrix of constants (n x 1) described in equation 


(4) for the system in Figure 3-1, 
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Figure 3-1. Flow Diagram of the system. 
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D =}. (4) 


B is an (1 x n) row vector, described in equation (5) for 
the system described in Figure 3-1. 

в- (Б, Б2...51 (5) 
The solution to (2) is given by: , 

1 ا 

ОБЕ Е EE? (6) 
Where 

e(t - ζω) = A Jm = ο, 


Solving the differential equation (2) for the discrete solution, 


X(k + 1) = @X(k) + DEL U(k) (7) 
Where 

to) (8)‏ 7 تپ 

DEL = MES t- t1) D at (9) 


Note that k(t - ا‎ is represented by К. 

Titus, in reference [10] developed a digital computer 
program called PHIDEL, This program provides a solution to 
equations (8) and (9) and will be used as a subroutine to 
obtain pand DEL. A listing of the PHIDEL program is provided 
in Appendix one. 

Plant Control 

It is desired to design a set of feedback coefficients to 

control the plant in accordance with a selected performance 


l U(t) is held constant over the interval (t - to) and is 
equal to. U(t,). 
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index. This set of coefficients will modify the plant states апа 
in effect will cause a movement of the Z-plane poles. This 
effect must be weighed against the feedback gains derived for 
the desired performance index. 

Let J(N) be defined as the performance index of a discrete 
sample-data system which will be minimized with respect to 


U(K), the control for the plant. 


N ë 
ΙΝ) Ξ Μία Σ [Xt(K) Q X(K) +r U2(K - 1)] (10) 
U(K) K = M 


Q is a (n x n) positive-definite symmetric constant matrix, 
and R is a positive scalar constant. 

The principle of optimality states that any portion of an 
optimal trajectory is also an optimal trajectory. Therefore 
equation (1) may be rewritten: 


J(N) = Min [X (N) Q X(N) +r U(N - 1) + J(N - 1)] Q1) 
U(N) 


If the gradient of J(N) is taken with respect to U(N - 1) and 
set equal to Zero, then an optimal control, U°(N - 1) may be 
determined. It is noted that J(N - 1) is not a function of 
U(N - 1). 


Using the above arguments, Titus developed the following 


algorithms: | 
At (k + 1) = -Át P(k) 4 
A P(k) A*r (12) 
Р(К) = yt (k)P(k - 1) Wk) +Q+ rA (kJA* (k) (13) 
Мк) = Ф+ ДАҢ(К) (14) 


Using the recursive relations of (12), (13) and (14), Titus [10] 
developed program OPTCON. This method is also discussed 
in references [3] and [4]by Titus, Strum and Demetry, and 
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in reference [6 ]by Ogden. A fortran listing of the program may 
be found in Appendix one. 
Digital Filter 

In the control discussion, the assumption was made that 
all the states of the dynamic system are observable states. This 
condition is certainly not always true and often the observable 
states may be measured only at the cost of some ambiguity due 
to the measurement noise, The inputs to the system such as 
the radar seeker discussed in Chapter Two will often be re- 
ceived ina noisy environment. 

Kalman in references [1] and [2] discussed these prob- 
lems and presented the theory for the desired filter, Schmidt 
5]; Titus, Demetry and Strum [4 ] and Jardine [7] have dis- 
cussed this problem and developed methods of implementing this 
problem on the digital computer. In an effort to maintain clafity, 
some of these developments will be presented. 

The digital filter will provide a best estimate of all the 
states by weighing the past information with the present obser- 
vable states. This weighting is performed with the knowledge 
of the environment (excitation) and the measurement noise, 
Although in many cases the noise is very difficult to describe, 
the assumption that white noise is present is in general true. 
Thus if white noise can be discriminated, the integrity of the 
measured signals will be improved. 

Filter Design 
For a single variable, the desnity function f(x) for a 


gaussian distribution is given by: 
2 
f(x) = 1 е-Ёх - u) 


(21) о σέ (15) 
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мћеге 


$) (ах = 2[9690] (16) 
where E implies expected value, 
and | 
бебдах = 1 (17) 


Assume that the systems measurement and excitation 
noise are white noise, where white noise has a gaussian dis- 
tribution and is spread uniformly over all frequencies. 

Noting that the input, U(k) to the plant is independent of 
U(k - 1), X(k) is independent of X(k - 1) and using the joint 
probability property of random independent samples, a density 


function for the state system may be written: 


F(X) = 1 QSAEX-XIPUOX-X] (8) 
ета’? рф ` 

мћеге Ў (к) is the best estimate of the states based on past 
information. 

Z(k - 1) = H X(k - 1) + V(k - 1) (19) 
Z is the noisy observable matrix, and V is the measurement 
noise vector BEE ру: 

Е [У] =0 (20) 

E [VVt ]- R (21) 
R is the covariance matrix of the measurement noise and P is 
the covariance matrix of the error. P is a symmetric matrix, 
(Pi; - Pj ). The diagonal terms of P are equal to the variance 
of each state (c ET and the off-diagonal terms are equal to 
correlation between the states, 

(E Û x, ° х,1- X; А x). 


кек es (22) 
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The expected value of the error is defined as the loss function: 


Λ А 
b= f x-S5tx- θεασά, αχ (23) 
^ 
Taking the gradiant of equation (23) with respect to X; 
VR L= fzxr(yz,3ax - 2X [к(х/2,Х)ах (24) 


Applying equations (16) and (17) and setting the result equal 
to zero 
^ 
VoL-2 E X/z,X - 2Х-0 


Thus: 

Х* = Max) = ELX/Z,X ] (25) 
X* is the best estimate of the states given the present noisy 
observable and the past prediction of the states. 
Filter Equations 

There are two methods of finding the recursive relations for 
the filter program. Method one assumes the random varables are 
gaussian and makes use of Bayes equation to derive an ex- 
pression for equation (25). Method two assumes a linear fil- 
ter, selects an algorithm for equation (25) and proves that this 
selection provides the best values for X*, 
Method One 
The covariance matrix of the observable error may be written: 

Рае ера (s); (25) 
Combining equations (19) and (26), 

P, = E[ (HX4+V - HX)(HX+V- 8011 
HEL(X- X(x- 3*3] + E [vvî ] 
Using equations (20), (21) and (22), 

P, = HPHt +R (27) 


l! 
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In order to find 57/60765 5162 Eeer Е 


applied to equation (25). 


A A 
X* = F(X/Z,X) = F(Z,X/X) F(X 
F(Z,$) (28) 


۸ 
Since Z(k) is independent of X(k), equation (28) may be 


written: 


A A 
F(X/Z ,X) = 2(2/Х) SC F (X) (29) 
F(Z) F(X) 
^ ^ ^ 
Since F(X/X) = F(X/X) F(X) , equation (29) becomes, 


F (X) 


F(X/Z,X) = F(Z/X) F(X/X) ei 


| F (Z) 
A 
Since X(k) is independent of X(k), (30) becomes, 


F(X/Z,X) = F(2/X) F(X) " 


F (Z) 
Following the form in equation (18), F(Z) and F(Z/X) are 
defined: 


3 ^ t =] А 
Е(2) = 1 е ЭК - 2) Р2 (2 - 2) (32) 
(21017 2]р,] ê 
Р(2/Х) = Е(У) = l еМ К” У (33) 


(2а)1/2 ірі 2 
Combining equations (18), (31) апа (33), 
Е(х/г Х) -де ?В 
Where 


A= _ (HPH + IS 
(2п0/ « |R| 2 |Ρ| 2 


= А ~ ۸ A E ^ 
p= vtR lv - x- Xtr!ix-3 «e (2 – 2) (НРИ! + в)71(2 - 2) 


апа 


Setting the gradiant of B with respect to X equal to zero. 
ot =] АЕ t -1 А 
VxB = М2(Х*=Х)°Р 7° - 2(2-2)(НРН +В) %7у(2-2) - 0 
біпсе 


А 
Vy (7-7) = V y (HX + V) = H, the algorithm, equation (35) 
28 


may be written. 


ο PR РН НМ ш 2) (35) 
Let the filter gain matrix G be defined. 

С = РН (HPH' + К)-! (36) 
Then the algorithm, equation (35) may be written: 

X* = X + G(Z - 2) (37) 
Method Two 


An alternate method of deriving the recursive relation for 
the gain matrix assumes a linear filter. The trace of the co- 
variance matrix is given by: 

L-E((x- 3x - £j (38) 
It will be shown that if L is minimized with respect to the gain, 
then the gain will be given by equation (36). 

Substituting equation (37) into equation (38), 
L= E[(X-X)(X-X)t] - ЕС -Х (2-Й 19 - | 
ECG(z-2) (K-X)t] + EL G(z-2) (z-2)*c* y (39) 
Since Z = HX + V 
Then Z2 E[ HX» V] - HX 
Assuming that E [ (x-X)Vt Jis equal to Zero, equation (39) may 
be written: 

L-E[QG-X)(x-3)t]- EC(-X) (c-3!] Híc* « GE(WWEj 

- GH E [(x-X%) (X-X)t] + GHE [ (x-X)(x-X)tutet (40) 
Substituting equations (21) and (22) into equation (40), 

L=P - pHtct - GHP + GHPHtct + Gret (41) 
Taking the gradiant of L with respect to G and solving for G. 

G = PH'(HPH't + R)7! 

Filter Programs 

Jardine, Titus, Demetry and Strum have designed programs 

to calculate the filter gains and the covariance matrix of error. 
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Such a program, (Program Filter) is listed in Appendix one. 
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CHAPTER IV 
DIGITAL SIMULATION 

If a digital filter can be obtained, then a method must be 
developed to evaluate this system, This chapter will describe 
the simulation of the missile tracking system, the digital pre- 
dictor, the control, and the kinematics. The results of this 
simulation for a particular missile will be discussed and tra- 
jectories from computer runs will be presented in Chapter Five. 
Missile Guidance 

Figure 4-] is a block diagram of the digital simulation. As 
described in previous chapters, the missile system tracks the 
target, obtaining an input for the missile guidance, This signal 
is proportional to the angular rate of change of the line-of-sight 
angle. In this simulation it is assumed that white noise is 
added to this signal. 

The details of the missile guidance simulation are des- 
cribed in Appendix Two. This portion of the simulation de- 
Scribes the tracking, control and dynamics of the missile as 
a transfer function relating the line-of-sight angular rate to the 
missile flight path angular rate. This transfer function is then 
formed into an F and D matrix as described in Chapter Three. 
Program Phidel is employed to obtain the тапа A matrixes for 
the state difference equation. 

The output of the missile guidance (y) is considered the 
observable for the predictor. Gaussian measurement noise 
with standard deviation (oy) is added to the observable. The 


output (y) is also an input to the kinematics. Figure 4-2 
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Figure 4-1. Block Diagram of Missile Simulation. 
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U(k) + | UNIT X(k) 
а с: а 


Figure 4-2. Discrete Flow Graph of the Plant. 
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describes the digital simulation of the guidance. The following 
fortran listing updates the discrete difference equation as dem- 
onstrated in Figure 4-2. 
С THIS SECTION UPDATES THE STATE VECTOR X 

CALL RNDEV(NUNIF, DEV) 

W=SIGW*DEV 

CALL PROD(PHI,X,KN,KN,1,TEMP1) 

DO 803 1=1,KN 
803 TEMP2(I,1) = W*DEL(I,1) 

CALL ADD(XS ,DINP,KN,1 ,TELP) 

CALL PROD(AT,TELP,1,KN,1,TELP1) 

CALL PROD(DEL, TELPL, KN, 1, 1, TETE?) 

CALL ADD(TEMP1 ‚ТЕМР2,КМ,1,Х) 

CALL ADD(X,TELP2 ,KN,1,X) 


Digital filter- predictor 


The digital filter is similar to the filter described in 
reference [4 ]. The gain matrix (G) is evaluated each sample 
instead of using the steady-state values of the gain. The dis- 
crete flow graph of the filter is shown in Figure 4-3. 

One input to the filter is the noisy observable. This is 
the şum of the angular rate of the missile flight path (у) апа 
the measurement noise. This noise is assumed gaussian with 
mean zero and variance 0,. This measurement on the missile 
could be made by use of rate gyros. 

The other input to the filter is the deterministic input. 
Here this would be the angular rate of the line-of-sight (с) 


and the best estimate of its derivatives. In the simulation 
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Figure 4-3. Discrete Flow Graph of the Filter/predictor 
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these are obtained from the kinematics. In the missile this 
signal would be measured by the tracking system. Note that the 
present value of the deterministic input DI(k) is used as the 
input to the plant while the past value of the input DI(k - 1) is 
used in the filter. The reasoning is that the filter is taking the 
present value of the observable Z(k) and therefore must use the 
past value of the deterministic input DI(k - 1) which caused the 
observable Z(k). 
The fortran statements for the filter are listed below. 
С THIS SECTION CALCULATES XS, THE BEST ESTIMATE 
C OF THE STATE VECTOR 

CALL PROD(H ,X,KP,KN,1,Y) 

DO 10 I=1,KP 

CALL RNDEV(NUNIF , DEV) 
10 V(I,1) = SIGV(I)*DEV 

CALL ADD(Y ,V,KP ,1 ,Z) 

CALL PRO DIERW xo, KN ,KN,1 ,TENMPI) 

DO 111=1,KN 

DO 11 J= 1,KN 
11 TEMP2(I,J) = -TEMP2(I,J) 

CALL PROD(TEMP2 ,TEMP1,KN,KN,1,TEMP3) 

CALL ADD(TEMEI TEMES EN 1P TEMES) 

CALL ADD(XS,DINP,KN,1,TELP) 

CALL PRODAT TELE IKN, 1 TEMPI) 

CALL PROD(DEL,TEMP1 ,KN,1,1 ,TELP) 

CALL PROD(TEMP2 ,TELP,KN ,KN,1,TELP1) 

CALL ADD(TELP 4 770 

CALL PROD(G Z. KN KP, 1 TELEP2) 
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CALL ADD(TEMP3,TELP1 ,KN,1 , XS) 
CALL ADD(XS , TELP2 ,KN, 1, XS) 
The Kinematics 

The kinematics combines the flight path angular rate of the 
missile, the missile speed (assumed constant) and the target's 
velocity vector to determine their effect on the line-of-sight 
angular rate (0) and the range rate (R) Е 

Figure 4-4 demonstrates the sign conventions used in this 
simulation. The trajectories which will be discussed later in 
Chapter Five use this same notation and sign convention. 

Figure 4-5 demonstrates in block diagram form the kine- 
matics. As noted in the figure, by-products of the kinematics 
are (1) the position of the missile and target, (2) the range and 
(3) the sign and magnitude of the range rate. 

By addition of a few fortran statements, adjusting the 
target's velocity vector, the target can be maneuvered as a 
function of the missile trajectory in range (evasion) or on a 
predetermined course (attack). 

The fortran statements for the kinematics are listed below. 
C THIS SECTION DETERMINES THE EFFECT OF THE 
C PLANT OF THE KINEMATICS 

GAMDOT = X(1,1) 

GAMMA = GAMMA+GAMDOT#DT 

XMDOT = VM*COSF (GAMMA) 

YMDOT - VM*SIN..F (GAMMA) 

YRDOT - YTDOT-YMDOT 

XRDOT - YTDOT-YMDOT 

RDOT = ((YRDOT*SINF (SIGMA) )+(XRDOT*COSF (SIGMA))) 

RDSIG = ((YRDOT*COSF (SIGMA))- (XRDOT*SINF (SIGMA))) 
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Figure 4-5, Block diagram of the Kinematics. 
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RLOS = RLOS+RDOT*DT 

DSIG = RDSIG/RLOS 

SIGMA = SIGMA+DSIG*DT 

XMZ = XMZ+XMDOT*DT 

YMZ = (۸2+7٣ 

XTZ = XTZ+XTDOT*DT 

YTZ = YTZ+YTDOT*DT 
The Control 

In this simulation a control (АТ) was selected through the 
use of program "OPCON" which minimizes the states X(k), ie, 
"bang-bang control". Other forms of control are also avail- 
able such as placing a limit of the fuel or energy expended, 
placing a limit on the amount of acceleration allowed, or plac- 
ing a limit on the magnitude of the angle and angular rate that 
can be measured. These options are recommended for future 
study and would, of course, be necessary in making a more 
complete study of the missile systems. 

For this simulation the weighting vector (At) becomes a 
constant vector which is imposed on the difference between 
the best prediction of the states X*(k), and the deterministic 
input, DINP(k) to determine the control to be applied. 

Figure 4-6 is a discrete flow graph of the entire simulation. 
Ihe fortran listing of the entire program may be found in Appendix 


One. 
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Figure 4-6. Discrete flow graph of missile simulation. 
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СНАРТЕК У 
SIMULATION RESULTS 

15 making a study such as the evaluation of a missile guid- 
ance, itis quite difficult to decide which parameters to vary and 
which to hold constant. It is even more difficult to decide which 
curves hold the most meaning as a measure of success. Since 
the ultimate goal is to hit the target, trajectories and confirming 
print-out were used as a primary measure. Once a hit was ob- 
tained with a set of noise variances, the target trajectory was 
varied to see the effect on the missile. The noise was also 
varied to study the amount of noise that could be tolerated. 

The measurement noise was found to be the most critical 
value. A standard deviation of measurement noise of 0,1 radians/ 
sec was found to be the upper limit. Above this value the 
missile had some control, but the missile trajectories were far 
from desirable. Operation near the limit of measurement noise 
caused large ambiguities of the target information at the be- 
ginning of the flight. The target information improved as the 
range decreased, however, and in most cases the missile was 
able to capture the target. 

Curves of the normal acceleration of the missile, the filter 
gains, the range, and the control versus time were considered 
and are shown in the results where these values appear to be 
critical or of interest. 

This simulation assumed that the missile propulsion has 
already burned out and the missile is traveling at a constant 
velocity (V,,). This simulation considers only the coplanar 


Situation, The program could be easily expanded to two channels 
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Figure 5-1 INITIAL CONDITION OF TARGET 


Target velocity direction: A,B,C,D 
Missile: Same initial position for all runs 
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and crosscoupling could be considered. The trajectories may be 
considered as yaw or pitch maneuvers, If they are considered 
pitch maneuvers, the effects of altitude variations and gravity 
must be considered. Otherwise it is assumed that the missile 
has the same dynamic reaction to a change in either plane. The 
gravity consideration was not considered in this simulation. 

The target was made to under-go several maneuvers: (1) а 
turn, (2) a straight line course, and (3) a left/right/left turn 
(evasive course). Figure 5-1 is a graphic representation of the 
missile and target initial conditions where A, B, C, D represents 
the target velocity vector for the four cases, The numerical values 
of all the parameters involved are listed in Table 5-1. The diff- 
erent maneuvers of the target are shown in Figure 5-2. A par- 
ticular computer run will be designated (A-1-2-0,1-0,1). The 
"A" implies the initial conditions of the target (position and 
velocity direction). The first "1" implies that the magnitude 
of the target's velocity is 1000 ft/sec. The "2" implies that 
the target is flying a straight line course (maneuver two). The 
last two numbers are the standard deviation of the excitation 
noise and the measurement noise respectively. Note that all 
angles in this program are given in radians. 

Figure 5-3 to Figure 5-26 are presented below to display 
the filtered missile capabilities, Listed on each figure are the 
range at which the program was terminated (R <100 ft.) and the 
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APPENDIX I 


THROUGHOUT THIS THESIS COMPUTER PROGRAMS WERE USED, MADE 
REFERENCE TO OR EXPLAINED. THIS APPENDIX WILL LIST- THESE PRO- 
GRAMS AND GIVE A BRIEF EXPLANATION OF THE DATA REQUIRED. 


PROGRAM FILTERe 


THIS PROGRAM CALCULATES THE OPTIMUM STEADY-STATE GAINS 
G, FOR THE DIGITAL FILTER. THE DATA CARDS ARE EXPLAINED IN 
THE INITIAL COMMENT CARDS OF THE PROGRAMe 


PROGRAM FILTER 

D1 ORDER OF SYSTEM IN I2 FORMAT 
D2 SAMPLING INTERVAL IN 8F104,0 FORMAT 
D3 F MATRIX BY ROWS IN 8F10e¢0 FORMAT 
D4 D MATRIX BY COLUMN IN 8F10.0 FORMAT 
05 VARIANCE OF EXCITATION NOISE SIGWSQ IN 8F10.6 FORMAT 
D6 R COVARIANCE MATRIX OF MEASUREMENT NOISE IN 8F10¢6 FORMAT 

PHI SYSTEM TRANSITION MATRIX 

DEL DISTRIBUTION MATRIX 

G OPTIMUM GAIN MATRIX 

H OBSERVABLE MATRIX 

P BEST ESTIMATE OF ERROR COVARIANCE MATRIX 

Q EXCITATION NOISE COVARIANCE MATRIX 
DIMENSION P(12912) 9Q(12912)9H(12912)9R(12912) 5G(12912) sPHIT(12912 


1sPHI(12s12) DEL (12)sDELDELT(12 12) F DELS(12+612)sDELST (1251 Dp 
2Х(12)»ХНАТ(12)»УНАТ(12) »ΡΝΕΝ(12»12) 


READ 33» N 

FORMAT (I2) 

READ 2,DT 

FORMAT (8F10-0) 

PRINT 1003 

FORMAT (1H1) 

CALL PHIDEL (PHI »DEL»NoDT) 

DO 1001 LL=194 

READ 20sSIGWSQ 

FORMAT(F10.,6) 

DO 1001 LK=194 
INIALIZE MATRICES FOR EACH RUN 
DO 10 I=loN 

DO 10 J=lsN 

PNEW(IsJ)20.40 


PC(IsJ)2040 
G(I9J)=0e0 
Пита) 20. 0 


DELDELT(I$J)20,0 
ВЕР ЕЁ: У} 2930 
DELST(I9sJ)204,0 


READ 2001,R(191) 70 





PROGRAM FILTER CONTINUED | 


2001 FORMAT(8F 10,6) : 
DO 30 108 Е 
30 02 5(1»1)=0Е (Т) d 
CALL TRANS(DELSSN,N» DELST) 
CALL PROD(DFLS,DELSTS.,N»N»NS,DELDFLT) 
DO 40 T=leN 
DO 40 J=leN 
40 Q(IsJ)5SIGWSQ*DELDELT(I»9J) 
SIGWSSQRTF(SIGWSQ) 
PRINT 10045.SIGW 
1004 FORMAT(/ 10X»5HSIGWs /5$E20,8) 
Р(191)-.02 
Р(102) =0,0 
Р(2+1)=0%0 
Р(2»%2):1.0 
PRINT 2002»R(1»12»((QUV I4) 5J21, ND s I1» ND) s (CPC Io J1) о Ј= ]ј о) 9» I2 15 N) 
2002 FORMAT(///20Xs8HR(1»1)* »F10646/20X.s8HQ( Ios J) 5 „4710.6/2О0Х „ВНР( |.) з 
1 ,4аҒ10.6) | 
Н(191):1. 
Н(1»2)0.0 
PRINT 700 
700 FORMAT (//5X93HG11 » 7Х + 38012 +77 2 , 38021) 7X9 3HG22 9 7X 4 3HP 119 7X» 3HP 129 7X» 
1 3НР21,7Х.,3НР22/) 
DO 1000 КК-1»40 
CALL GP(HePHI sPeQeRe20l sGePNEW) 
PRINT ВОО((0О(1+ Ј)• Ја је М) о] =1]1•М ) • ((РМЕМ (1 • Ј) eo J=l ON) oT El ON) 
800 FORMAT (10F106%5) 
DO 11 I21»N 
DO 11 J=1 +N 
11 P(IsJ)sPNEW(I9J) 
1000 CONTINUE ` 
PRINT 1005 
1005 ҒОВМАТ(ІН1) 
1001 CONTINUE 
END 
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SUBROUTINE GP(HePHI sPeQeReoKNoKP oGePNEW) | 
ОІМЕМ5ІОМ Н(12912)»РН1(12912)»Р(12912)»0(12912)9Б(12912)»6( 12912)». 
1ІРМЕМ(12912)ҺНТ( 12,12) о ТМ1 (12912) 5» TV2(12912) 
CALL TRANS CHeKP eKNoHT) 
CALL PROD(P yHTy yKN | KN KP TV1) 
CALL PROD( (Hy yTVl1e KP КМ „КР ,ТУ2 ) 
CALL ADD(TV2 FR KP KP  TV1) 
CALL RECIP(LKP,,00000000000019»TV1, TV2 »KER) 
IF(KER-2) 101»1105101 

110 PRINT 111 

111 FORMAT (5HKER=2) | 

101 CALL PROD(HT ə TV2 | KNe KP | KP e6fTV1) > 
CALL РКОО(Р»ТУ,ҺЭКМЭКМӘКР» 0) š 
EALL PRODI Har KB KN KN TV) : 
CALL PROD(GoeTV1 oKNoKPoKNo TV2) 3 
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PROGRAM FILTER CONTINUED 


DO 102 [=1›КМ 

ا۳۷ 1ل 102 ОО‏ 

TV2(I sJ)=-TV2 (I J) 

CALL ADD(PsTV2sKNsKNoTV1) 

CALL PROD(PHIsTV1sKNsKNoKNoTV2) 
CALL ТВАМЅ5(РНІ, КМ, КМ» ТМ] ) 

CALL PROD(TV25TVIS.KNsSKNsSKNSPNEW) 
CALL ADD(PNEW 9Q5KN sKNsPNEW) 

END 


SUBROUTINE TRANS(AsN9Mo9C) 
DIMENSION A(12912)5sC(12912) 
ОО 153 [ = 1›М 

DO 153 ром 

С») = А(1•Ј) 

END 


СВВОУТІМЕ РЕСІР(МЭЕР»А»Х>КЕВ) 
DIMENSION А(12»12)»Х(12»12) 
ОО 1 Із1»М 


ОО 1 Js1lsN 

Х(1»Ј)20%0 

DO 2 Kz1»N 

Х(К)»К)< 1.0 

DO 34 Lz1,N 
КР=0 

2-00 


DO 12 K=LoN 

IF (Z-ABSF (A(K9L)))11912912 
2-АВ5Ғ(А(К»()) 

КР = К 

CONTINUE 

IF(L—KP)13 20520 

DO 1& J=L+sN 

Z=A(L +J) 

A(LsJ)=A(KP9J) 

А(КР+Ј)=2 

ОО 15 18ل‎ 

2-Х(10».) 
Х(1»9.2)-Х(КР,,) 

Х(КР+Ј)=4 
ТЕ(АВБЗЕ(А( οἱ) -ΕΡ) 50.50.20 
ІЕ((-М)231»34»34 

LP1=L+1 

DO 36 К-іІР1»М 

ІР(А(К»()) 32»36,32 
RATIO=A(K+L)/A(L +L) 

DO 33 J=LP1sN 

А(К».)) <А (Ко .)) -ВАТТОХЖА (| ».)) 
ба 38 Ј=1»4 
Х(К».4/)-Х(К».)-КАТІО%Х(1»>.,) 
CONTINUE 
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43 


50 


1003 
399 


3991 


400 


500 


401 


PROGRAM FILTER CONTINUED 


CONTINUE 

ОО 43 1=1,М 
II=N+1-I 

DO 43 J=loN 

6-0.0 

ΙΕ(ΙΙ-Ν)41 43943 
11Ρ1Ξ1 15411 

DO 42 ΚΞΙΙΡ]1.Ν 
5=25+А(11К)%Х (Кој) 
Х(11»/3<(Х(11»./)-5)/А(11,11) 
КЕК=1 

RETURN 

KER=2 

END 


SUBROUTINE PHIDEL( PHI sDELeNoDT) 
VALID ONLY FOR A CONSTANT F MATRIX 
DIMENSION F(12912)ePHI (12912) ,TERM( 12912) sWORM(12912) 


1» DEL(12) sDELM( 12012) sTELM(12.12) sDELP( 12012) 20(12)} 


КЕАО1, ( (Е (ТК,„ТС) ,ТС- |» МН) о 1 К= 1,03) 

ЕҒОВМАТ ((8Ғ100)) 

READ 1 (0(1)»іІг1» М) 

РЕТМТ 399ЭОТ»(ІК(ІР, ІС)» ІСг19М)» ІК1» М) 
ЕОКМАТ(///З3НОТ:»)1Е543///УЭ7НҒ(1,42/3х/,)((8Ғ82))) 
РЕТМТ 3991 (0(1)»1=1•%М) 

РОВМАТ ( ///75Н0(1)=/(82Е8+2)) 

ΝΕΙΝΑΙ 51 

ІМ-0.0 

DO 400 IR=1,N 

DO 800 ІС-1,М 

ТЕКМ (1К›][С) =ОеО 

МОЕМ(ІР»ІС)-0.0 

ТЕКМ (1В›]К) =] еО 

ТЕ! М( ТК» | С) -ТЕКМ( ТК, i C) *DT 

ПЕР ТК» 1 С) -ТЕ! М (ТК, ТС) 

ОЕСМ (ТЕ ТС) 100 

DEL(IR)30,0 

РНІ (ІР, ІСҘ)гТЕКМ(ІР, ІС) 

ТМ= ] е0+ТМ 

00 500 1® =] ЭМ 

00 500 ІС-1»М 

ОО 500 JN=1 +N ۱ 
ОЕСМ(ТК.,ТС) -ОЕ! М (ТК, ТС) -ТЕ: М (ТЕ, „М) ТЕ (ЈМ, 10) #07 /(ТМ+] 60) 
WORMUIR9$IC)ZTERMCIR9$ JN) *FCJNS IC) XDT/TM*WORM(IRS9 IC) 
DO 801 ІКг1,М 

DO 401 IC=leN 

ТЕКМ (ТК, ТС) + ИОКМ ( ТК, ТС) 

ТЕ СМ (ТК9.ТС) = ОПЕ М(1К+ 190) 
DELP(IRSIC)zDELP((IR9$? IC) * TELM(TR$IC) 
пете БЕІН Т пете ТЕРМ ЕРУ Є) 
DELM(IRsIC)s0,0 

WORM(IRsIC)s0,0 


75 


600 


503 


152 


191 


PROGRAM FILTER CONTINUED 


АВС-0.0 

ОО 21:1»М 

DO 2. Ξ1«Ν 

AA=TERM( I J) 

AB=ABSF (AA) 

IF (ABC~AB) 359392 

ABC=AB 

CONTINUE 

IF(0.000000005-ABC)5 556 

GO TO 4 

PRINT 502s((PHI(IRsIC) » 1221 М ) , Т2=1»+М) 
РОКМАТ ( / 1 11/9Х + ВНРНІ( 1973) /// (6Е208 ) ) 
DO 600 1=1, М 

ОО 600 К-1»М 

DO 600 Jz1»N 
DEL(I)ZsDEL(CI)TPHI(CIS J) *XDELP( JS, K)*D(K) 
PRINT 503 (DEL(I)sI=19N) 
FORMAT(////9X s6HDEL(1)//(6E20.28)//) 
END | 


SUBROUTINE ADD (А»В»М»М»С) 

DIMENSION A(12912)sB(12912)sC(12s12) 
DO 152 Іг1,М 

DO 152 10ل‎ 

ССТеју = АС Рој + B(I J) 

ЕМО 


SUBROUTINE PROD (AsBoNoMolsC) 
DIMENSION A(12912) 9B(12912)9C(12912) 
DO 151 1-1٣0 

DO 151 J=1sL 


С(1•Ј) 29, 

DO 151 K = 19M 

C(IsJ) = C(IsJ) + A(1sK)#%@B(K sJ) 
ЕМО 

ЕМО 


РКОСКАМ ОРТСОМ 
DATA CARDS 


D1 CASE THE COMPUTER IS TO EXAMINE. 
CASE 1 Q = Is R = O 9 Q* = 0 
CASE 2 Q = Os R 1 9 Q* = 0 
CASE 3 Q = 1, К 1» Q* EMPLOYED 
О2 М = ORDER OF SYSTEMs NSTAGE = NO, OF ITERATIONS 
D3 Q-WEIGHTING MATRIX 
D& F MATRIX 
0 6 +۳۴ 
D6 DT~SAMPLE PERIOD 


7% 


о Су СУ СУ СУ СУ CV CV CY CY су су 


PROGRAM OPTCON 

DIMENSION PHI (12012) sPSIT (12012) sPl(l2s12) sDEL(12) eAT(12)¢ 
10М(12»+12)»+»0(12+12)»+»Х( 900)» 17172 (12), БЕМ(16), ЕМ( 12) ,1НМ(12+ 12) 
1 Ү1(900)»Ү2(900)»Ү93(900) 


THIS PROGRAM UTILIZES A COST FUNCTION, J(N)=MINIMUM(SUM X(N)T*Q%*X(N)+ 
SUM R*#U(N-1)**#2), AN UNLIMITED NUMBER OF ITERATIONS MAY BE MADE AT 

A COMPUTATION RATE OF 2000 PER MINUTE AFTER THE PROGRAM HAS BEEN 
COMPILEDe THE OUTPUT OF THIS PROGRAM IS THE FEEDBACK GAIN MATRIX» 

A TRANSPOSE». THE FOLLOWING RECURSIVE EQUATIONS WERE DERIVED USING 
DYNAMIC PROGRAMMING» 


AT(K)=~-(DELT*P(K-1)*PHI)/(DELT*®P ( K-1]1) *DEL+R) (1) 
PSI(K)zPHI*DEL*AT(K) (2) 
P(K)sPSIT(K) *PCK-1) *PSI (K) «Qe RYA(KO)*AT(K) (3) 


Р(0)0»АТ(0)-0» Р51(0)-0 
EQUATIONS 1» 2» AND 3 CONSTITUTE THIS PROGRAM 


CALL IN DATA AND INITIALIZE 
DO 1111 1=1,3 
READ 30sKASE 

30 FORMAT(I1) 
READ 1 sNeNSTAGE 

] FORMAT (8110) 
READ 2»К»0Т 

2 FORMAT((4F2040)) 
READ 2 ((Q( I oJ) sJ=1loN) э1 =] эМ ) 
PRINT. THE DATA READ ІМ» 
PRINT 100 


100 FORMAT(1H1) 


PRINT 32, „КАЗЕ 
32 ҒОЕВМАТ(9Х%5НКА5Ес 213) 
PRINT 3eNeNSTAGE 
3 FORMAT (//9Xs2HN=013020X0 7HNSTAGE=013) 
PRINT 4eReDT 
4 FORMAT (/9X o2HR=0F15e1192Xs 3HDT=5F 15011) 
PRINT 9ef (QI 9J) sJ=loN) sIT=l oN) 
9 ЕОКМАТ(/9Х »7Н0(1+ Ј)=2/(2215•11)) 
CALL PHIDEL(PHIsDEL,NsOT) 
ОО 5 15.1.Ν 
DO 5 Jz1,N 
GM(IsJ)20,0 
EM(I)20,0 
ЕМ(І)-0.0 
5 Р(1+Ј)=0(1+Ј)) 
INITIALIZE P(O) FOR CASE TWO 
IF(KASE-2) 355936035 
36 Р(1+1)=1%0 
Р(2,2):1.0 
Р(3,3):1.0 
35 CONTINUE 
PRINT 19 
19 FORMAT (//2X s6HNSTAGE »› 4Х. 7НАТ ( ] ›]1 ) »ЗХ», 7 НАТ (12 ) »4Х»6НР( 1), |) » 
}&Х»@НР (1,2) ›,4Хоб6НР (2 3] } »4Хоб6НР (2.2 ) ) 
CALCULATE AT(J) 
DO 22 KK=lsNSTAGE 
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PROGRAM OPTCON CONTINUED 


ОЕМ-0.0 
006 1+1,М 
006 Jz1»N 
6 ‘EM(I)=EM(1)4+DEL (J) *#P( 451) 
008 Іг1»эМ 
DOT J=loN 
7 FM(I)=FM(1I)+EM(J)*PHI(J51) 
8 DEN=DEN+EM( 1) *DEL(I) 
DEN2-DEN-R 
0010І-1»М 
AT(I)=FM(1)/DEN 
ЕМ(І)-0.0 
10 ЕМ(І)-0.0 
CALCULATE Р5І (К»1І».) 
0013І1»М 
0013./=1 эм “ 
13 PSI(IsJ)=PHI(1sJ)+DEL(1)*AT(J) 
CALCULATE P(K»IsJ) 
ОО 15 [=1›М 
DO 15 J=loN 
DO 15 L=loN 
15 GM(IsJ)=GM( IT sJ)+PSIT (Lo TI) *P (LJ) 
DO 17Іг1»М 
DO 17/:1»М 
DO 16 L=loN 
16 HM(I+J)=HM(I+J)+GM(I +L )%PSI(L , J) 


CASE 1 TERMINAL CONTROL» OMIT Q(I»J) IN ЕОЈЏАТТОМ РОК Р(1»Ј) 
CASE2 MINIMUM FUEL OMIT Q(IsJ) IN EQUATION FOR Р(1+Ј) 
IF(KASE-2) 31931533 | 

31 Р(1»Ј)=НМ (1) Ј)+ВХАТ( 1) ЖАТ(Ј) 
GO TO 17 


CASE THREE MINIMIZATION OF TIME AND FUEL 
33 P(IsJ)=HM(I+J)+Q(I+J)+R*AT(I)*AT(J) 
17 HM(IsJ)20,0 
0018 1.=1›М 
0018 )-1,М 
18 GM(IsJ)=0-0 
PRINT 20sKKsAT(1) sAT(2)9P( 191) 9P(192)sP( 2591) sPl292) 
20 FORMAT (159(6F10¢6) ) 
22 CONTINUE 
1111 CONTINUE 
END 
SUBROUTINE PHIDEL( PHI »DEL»sNoDT) 
VALID ONLY FOR A CONSTANT F MATRIX 
DIMENSION F(12912)sPHI(12912) 5TERM( 12912) sWORM(12912) 
ls DEL(12) sDELM( 12912) sTELM(125912) »DELP(12912) ,D0(12) 
1 FORMAT ((8F10-0)) 
READLIs ( (FC IRsIC) sIC=l9N)sITR=l15N) 


READ 1 (0(1)» із1,М) 
1003 PRINT 399ЭОТ»((Ғ(ІР»ІС)» ІС-1»9М)»,ІК-1»М) 


399 ҒЕОКМАТ(///3Н0Т-91Ғ5.3/ ЭТНЕ(1».2/)-/»((2Ғ8.2))) 
78 


3991 


. 400 


500 


401 


600 


12 


PROGRAM OPTCON CONTINUED 


PRINT 3991 (DCT )eT=10N) 

FORMAT (/ 5НО (1) =/(2ЁЕ8е2 ) ) 
МЕТМАСЦ = ] 

۲00 

DO 400 IR=1,N 

DO 400 ICz1,N 

ТЕКМ ( ТК. ! С) 20 .0 

НОВМ ( ТЕ. ,ТС) +0.0 

TERM( IRs IR) =120 

TELM( IRs IC) HTERM( IRS IC) ОТ 
DELP(IRsIC)=TELM( IRS IC) 

ОЕ!: М (ТВ. ,„ТС) 20 .0 

ОЕ (ІР) 0.0 

РН! (ТК „ТС ) <ТЕКМ (ТК „ТС ) 

TM=1e0+TM 

DO 500 IR=1 +N к 
ОО 500 IC=1 +N 

DO 500 JN=1 N 
DELMCIRSIC)ZDELM(IR$S IC)-TELMCIRS JN) *FCJNS IC) *DT/( TM*160) 
ИОВМ (ТВ ТС ) -ТЕКМ (ТВ, АМ) КЕ (МТС) “ОТ / ТМ + МОВМ ( ТК. | С) 
DO 401 IRz1,N 

DO 401 ICs1l»N 
TERM(IR$IC)ZWORM(IRs$ IC) 
TELM(IRyIC)=DELM(IR y IC) 
DELPCIRsIC)I=DELP(IR»sIC)+TELM(IRs IC) 
PHICIRs IC)=PHI(C IRs IC)+TERM(IR,IC) 
DELM(IRsIC)=000 | 
WORM(IR9»1C)20.,40 

АВС-0.О 

ОО 211»М 

DOr 2J=1 oN 

АА=ТЕКЮМ ( 15.) 

AB=ABSF (AA) 

IF ( ABC-AB) 3,3,2 

АВС=АВ 

CONTINUE 

ІЕ(0.000000005-АВС)5»%5»6 

GO TO 4 

РБІМТ 11» ((РНІ(1»24/)»)21»М),1Іг1»М) 
РОВМАТ ( /ЭХ» ЭНРНТ ( 1») =/(2Е15.11)) 
DO 600 I=leN 

DO 600 Кг1»М 

DO 600 Ј=1, М 
DEL((I)SDELUI)*PHI(CIS J) XDELP( JS K)*D(K) 
PRINT 12e(DEL( I) oI =lsN) 

FORMAT (9Xe7HDEL( I) =/(1F150e11)) 

END 

END 
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PROGRAM ESILE 
ПАО ВИНО АЛИ SIMULATE 
THE KALMAN FILTER AN 


2 0 
THIS PROGRAM 


ТЕ СИРЕ ОБЕ А CLOSES om 
یت‎ ۳۳0۴ ۱5۶۱۱۱۱۹۸55006 
۶ی۶۶"‎ ۹ ON 7۶ 
ПЕ ЕЕ ЕТЕК САТЕ т 
۶ ۷ R K 


lide PROGRAM SOLVES Mi 
КОРЕК (K 
Εκ πο κι ο 

ο ο πι D Gn ۰۶ 


Eege 


X(K* 1) 5PHI*X(K) DEL ( 
IS MEASUREMENT 


TERE № 
DINP IS THE 25 Рале 


DATA CARDS. 


А ЗЕЕ ИН 
Ва САМА = “FOCH 
ие 
EE Ee 
0516 = LCS ANG 
D3 RLOS = RANGE 
КЛ а ο | Re Oar 
ο ο [STR X 
“OOP Pp hb hia 
И 
ΥΤΡΟΤ = INITIAL 
04 INITIAL CONDITI 
D5 AMZ INITIAL X-C 
ОЕ ИЕ ЕЕЕ 
۰۰۰۰۰٣ 
РОВЕР ЕТ: = 
ОГ Е МАТК:Х 


DoD VECTOR 
Bo GRAPH E 8 
ЮО АТ ες ος S 
bs sles 
012 $16м = S 
519۰۷ = Si 


DIMENSION 
ШЕ ото) наа А 
ο ος 2.512) Τεν. 
= 12)» 01 МР (4 2: 2 зе 
4эХМ (900) ›УМ (900 } 5Р (12 


۱٣٤۷7٦‏ و 
) 


— ——— es 


USES VARTASLE 


ee эу... 


5 THE MISSTLE CONTROL sDYNAMICS » AND ΧΙ. 
D Sa, CONTROL IS APPLIED TO THE СОМ mim 


FILTER GAINS | 


НЕ LOOP ON THE OPTIMUM OB : ا‎ 2 ٣۷ ء٣‎ 
THAT STABLE CONTROLLER GAIN MATRIX НАЗ 
GASIS OF DESIRED RESPONSE ΑΝΘ Ελα... 

X EACH БАМРСЕНОМ ти > 

OPERTIES OF THE ANTICIPATED س٤۹‎ 






E FOLLOWING 0 ОМО 


GH )2DEL*AT(XS(K-1)—-DINPeK=- 1 ΠΠ 
0د ۲ ۶ؤ‎ P 
IS THE RANDOM ٣151 ۸۵8 


А ГУ ( 
4 


1 RANGE FOR HET 

IND INA TE VOR. TA. Gel 

FLOCTTY ІМ х ОТКЕСТТОКМ РОК ТЕ Паси 
OSSDINAYSEOGR TARGET 

IN У ОТБЕСТТОМОХ PHE TOR m 
М ΤΗΕ STATE УЧЕСТЬ 

ORDINATE JO EMI SSITLE 

ORV LNA T ο 6 


0 
wae 


cm 


/ 
مسد 


En 0 EN 


/ 


0 ۹٣ 


(ay 


1 
r- 
8 
e 
к=з 
zd 


— 


ооо О 


N FOR CONTE ) ۹۷ 
МЕ TWO 
D DEVIATION 
DENTATION 


OF EXCIUNM T Ор от 8 
OF MEASUREMENT NOISE 


12) 5 ІСУ(12)>9Ү(12,12)>2(12>12у)эРНГС ИИ 
312) 0(12»12) „ТЕМР 1 (12, 12) „ТЕМР2( 12 
TELP (12212) ТЕРА 12) 12) „ТЕ: Р2 ( 1 2» И 
) $ s АВР АДИС 
EN WED 


> 


80 


Су Фу 


MOO 


N W 


кз 


+ 
x 


O» 
KH 
© 
τ 


200 


131 


Ша 


5909 


PROGRAM MISSILE CERD 


ο ο о DELS 
ο ορ ο ТЕРҮ ООО 
BNET SOME ο... 


READ 


2,М,0Т 


БЕАО 1, ОАММА 5! ОСМА) АМОС, 99510 
р ۲ 0001 ۴د‎ ٣624٦ 0 07 

ری رو را ا ٹن e‏ 

READ 1 59KMZ5YMZ5VM 

FORMAT (GAR? 

о 

ΠΤ ο аи 

КЕЗО СР го МЕДА) 

тая Eu nang ieee Py 

EE ADS IJ 92-126! 

READ РУСАТ Тотев, 

ο ος Са Тре 2, 

ΕΛ По) 

PRINT 144 

BORMATVGASGMAM: GCS i CEG oAs CHAISRGE I sSAsSnYMISSILE sexs SRYIARGET 
DO 6000 = тум I 
DELCISI1)SDELS(J) 

ШІСІ” ТКАК6 СОРТ Ев] 

CALL РКОрСРЕ 02 те МЕ ОГ ) 

MO 200 Jy 

рО 200 ٣إ‎ 

О: ба 
и. 

Duos ud cA) 

ПИ о ЕМ ORDER Ms mscwesiwNUNBEN OF ὈΒ5ΕΠΝΗΒΕΕΟν 
КМ =N 

КР=1 

и 

Ее АЕ ТНАТ Νο DETERMINISTIC 
РОТ РЕА ЫР ОУ TO TIME = ZERO 

ОО 1251 іс1%кМ 

(5. Arte EI ۱ 
NON F= 1220703225 

т С LONS 

GALL PROD LHs As Kos Ried 79) 

Dio NN Ке яу 

Cee RNDEVLNUNTESDEZV) 

ЕТЕМ ОТ ОЕМ 

САП: ADO IYO V o KP s ze) 

ο ο И Пт) 

ΚΚΞῚ 

I TU 

KK1=0 

о qus οκ MATRIX 


| بت ہے ہے‎ s Ge =: ٔ ,ة6‎ ! 
CA = = Ех(А Хх Рег а 02 Ка НОР, а) 


Па то Пи СДС ЛА ES KS tHe BEST ESTIMATE 


) 


СУ СУ СУ 


MOD 


WW о» 


F2 r2 


ΕΠΗ MIT SSIAE CON NE 


ri) 


D 


Oe ae ΤΑ УЕ ٦ 

CALL PROD(H»XoKP5KNo19Y) 

ОО 10 1=1.КР 

САБО ККОЕ М МОЕ О 

И ЕО ОЕ 

СА АБО (У, УКР. 1.2) 

САКЕ РВОО( РИХЫ ШКУ o 5 1) 
САГЕ РКОО(О,Н КН „КР КА  ТЕМР2 ) 

DO 11 РЕШЕ 

ВОД реда с 

ТЕМРӘ(Т» І--ТЕМРОСТАРН 

CARU PROD(TEMP2 0 EMP 3) 
CALL ADUCTEMP IST PES TS 
CALL ADDS LIND SS 
CALE PROD (ATS I o PKN aaa 
СА РБЕОО(РЕ о 1 2МР1, СМЕ FELS) 
СА! РЕОО(ТЕМР2, 1ЕГР,КМ,„КМ» ТЕ Р1) 
САВВА С РЕБРО о а 

CAL L PROD ( Спа смо Ро 5 52) 
САША ТЕХ IPAPITELPIGENSIOXS) 

CALE ADDI ASS TELE 2s CA el sks} 


DINEM а Бо) بب‎ ο 
(1۰۰۷۷۲١2 و‎ 46 
ео пе 
ТМЕ (КК ) =Т 


КОЕ ЕРЕ Т О М АР А ۰۰ ۸ ۶ ۹۶۸٦ 


CALL RNDEV(NUNIF ,DEV} 

WzSIGW*DEV 

CALL PROD( PHI Xy KNS ANAM  ТЕМР 1) 
с 

ТЕМР2 ( ο ο ο... 

CALL ۸0001: ο LP) 

CALL РКОО(АТЬТЕЕЗ EN: з ИЕ CE 
CADET PROD (DELS Feree alale) 
CARLE E KN 1 X 


CALL ADD(XsT! 


THES Ες ΤΟΝ ΡΕ ΤΕΙ ос ПЕЕ СЕРЕС ТИ Со и е И 
ON THE KINEMATICS | 
Ок516-0516 
ОТЕК(ОКК)-ХА( 29» 
GAMDOT=X( isi) 
GDOT(KK)2GAMDOT 
GAMMA =GAMMarGA4OO TROT 
AMDOT=VM*COSF (GAMMA } 
YMDOT=VM*SINF (GAMMA) 
ХЕРО РО Те eer 
XRDO SATO ЕВ 
КОСЕ ОРО SES 
( 


GMAT 7)0 6 EE 
05 OE S 


GMÀ))-OXRDOTESINF(SIGMA) )) 
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de 
< 
A 
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RB = 
το. 


De 
О 
e 


Lë 
D 


ҥе 


ον ο Os 
КЕМЕСІ 


UD QD r2 rS a. d 
O U U F3 N Ον 


IS 


PROGRAM MISS! 


x<-<x <x < Dë 


: 
ιν 
м. 


ОК = 205 
PRINT 140 ΚΚ 
PRINT 14 یئ‎ ) ç Á 


REOS ВЕО 5 ρο τω 

Ρο ο Кро Ра мн СЪ 

О = Бл ОСТ б 1 
Pos. O=(D51G-0.316770 5 
MZZXMZtXMDOT*OT 

۷ = ٣/۱۸2 ۷۷۲۵م‎ ۲۴7 

N رر‎ 


АА 
PORMAT АУ 9) 
ҒОБВМАТ(2Х;6Е15.7; 
акта 

| 
ИНЕ ОСЕ oS. К 
PORMA CO USS ST 
ο επ οτι δες م۹‎ Э 
FINT bus 07 

EE EE 

3.6 8 0ے و و خ1 ر000 س16 
MARCOSA CuI OU ὅμως.‏ 
όσο‏ .۳ط 

E πο τος 
FORMA. De КООТ=,1=16з7 
РОКА ιο NLOSSYE.-O-T 
дедо и πω 607۳001. 
о орон 4. 
мм s 

РОБ (Те ТМ. ЗЕЕ. СЕ 
رح ہمت ا‎ 
ENS Ie 

LAUD dl 

e. срам AM YN aT Cv 
VC 3 

LS AUD i 

CAL. ORAWUSK A T Y i 9 V C 9 
ο) 

ОСОО МЕ РАО ри ром 
pus SE O> 
ËO о. =» 

DOSE) Је о 

С.) = Оз 

м 

ος νο στο >‏ وا ےتک 
END‏ 

м к Ао» гу» 
l 222227631. 
DOTI? TEL y 


O 


ا 


ΤΟΙ ΘΕΣ 


УМК}, 


Сә {У 


к” 


t— 


АУТА SIGMAS 


1190>999>0>0»0 


ү 


12167) 


9795909LAST } 


О»0»0+7»8»0»А57) 


CY CV CY CY CY O 


PROGRAM MISSILE CONTINUED 


ЈЕНС Те ΕΛΙΤ. το. 
ЕМО 


SUBROUTINE FILTER(Ns KP I PHI,sQ R +H P G) 
PHI SYSTEM TRANSITION MATRIX 
DEL DISTRIBUTION MATRIX 
G OPTIMUM GAIN MATRIX 
H OBSERVABLE MATRIX 
P BEST ESTIMATE OF ERROR COVARIANCE MATRIX 
Q EXCITATION NOISE COVARIANCE MATRIX 
DIMENSION P(12912) 9Q(12912)9H( 12912) »R( 12912) 56(12512)5PHIT( lage 
L»sPHI (12°12) sDEL(12),DELDELT(12,12) sDELS (12912) БЕСТ ٦7٣ 
2PNEW( 12912) 
CALL GP(H»PHI9»P,Q,R; N, KP, G,PNEW) 
ОО 11 Із1»М 
ΡΟ 11 .2/:1»М 
11 P(IsJ)=PNEW(I9J) у 
| ЕМО 





| 
5ОВКООТІМЕ ОР(Н»ЭРНІ»Р»С»КЭКМЭКР»ОЭРМЕМ) 


DIMENSION Н(12912)ЭРНІ(12912)»Р(12»912)»0(12»12)»К(12%912)»0(12>12 
IPNEW( 12912) sHT (12312) sTV1(12912)sTV2(12912) 
CALL TRANS(H»KPsSKNesHT) 
СА. РВОО(РАНТОКМ КМ ОКР» ТМ] ) 
CALL РВОО(Н»ТУІКР»КМЭКР»ТУ2) 
CALL ADD(TV2sRs KP, KP o TV1) 
CALL RECIP(KP$,,0000000000001*TV15» TV2 s. KER) 
IF(KER-2) 10151105101 
MO PRINT 111 
111 FORMAT(5HKER=2) 
101 CALL PROD(HTsTV2S.KNS,KP,.KP,TV1) 
CALL PROD(P,TVIlS.KNSKNSKP,G6) 
CALL PROD(HsPsKPsKNoKNoTV1) 
CALL PROD(GsTV1sKNsKPoKNo TV2) 
DO 102 Iz15,KN 
DO 102 او ےل‎ 
102 Ту2([›,.))=-Т\/2([›.)) 
CALL ADD(P,.TV2S.SKNsSKN»TV1) 
CALL PROD( PHI sTV15KNoKNsKNoTV2) 
CALL TRANS( PHI sKNsKNoTV1) 
CALL PRODI”(TV2sTV1, KN KN KN PNEW) 
CALL ADD(PNEW, IQ | KN KN + PNEW) 
END 


SUBROUTINE ТКАМ5(А»М»М»С) 
DIMENSION A(12s12)sC(12912) 


DO 153 |: = 1›М 
рО 152 2-1»М 

0383 Gri) = Ai indi 
END 
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34 
40 
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43 


PROGRAM MISSILE CONTINUED 


SUBROUTINE TRANS(AsNoMoC} 
DIMENSION A(12912) sCl1l2s12) 
DO 153 I = l»N 

рО 153 Ј=1,м 

C(J» I1) = А([»ј) 

ЕМО 


SUBROUTINE БКЕСІР(МЭЕР,А»Х» КЕН) 
DIMENSION A(12»12)5X(129»12) 
DO 1 T=10N 
DO 1 12ل‎ 
X(19J)=020 
ОО 2 К-1»9М 
Х(К,К)=10 
DO 34 Lz1»N 
КР=0 
72-060 
DO 12 K=L yN 
[= (2—АВБУР(А(К» ())) 11» 12» 12 
2= АВЅҒ(А(К, 1) ) 
КР = К 
CONTINUE 
ТЕ (| -КР) 1320 »20 
DO 14 J=L +N 
Z*A(L»J) 
A(L» 2) SACKP,J) 
А(КР•Ј )=2 
DO 15 Js1,N 
2=2Х( +J) 
X(LoJ) =X(KP oJ) 
Х(КРо.] } =2 
ІЕ(АВЅҒЕ(А( 51) ) ЕР ) 50» 50» 30 
ТЕ (| -М)31 34» 34 
ІРігі-1 
ΡΟ 36 К-іР1»М 
ЈР(А(К» 0) )32»36+ 32 
КАТІОзА(К»()/А(094) 
ОО 33 ӘзіРі»М 
А(К»Ј)гА(К» Ј)-КАТТОХА(  •Ј) 
DO 35 Jz1l,N 
XOK& SJ) ZXOK 9 J) -RATIOFX(L9J) 
CONTINUE 
CONTINUE 
DO 43 I=] sN 
II=N+1-I 
DO 43 J=1l sN 
$=0.0 
12 (11—4)8 1» 53» 83 
ІІР1=11+1 
DO 42 Кг-ТІРІ1»М 
SSS*AUDL TS KO * X CIL v d) 
XCIISJ)s2(XCIISJ)-SO)/ACI I SI I) 
KERz1 
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400 


БОО 
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PROGRAM MISSILE CONTINUED 


RETURN 

KER=2 

END 

УЈЗУВООТ ЊЕ РАТОВЕ НОРА У ВЕР bare 

DIMENSION (1212) РНТ 12» 12), ТЕКМ 12, 12), МОКМЕ ٦ 


Le 0ЕЁ( 12} .0ЕЁБМ( 12,12)» ТЕЁМ (12,12), РЕБРО 2512) PD D 


VALID ONLY FOR ۰ھ‎ 0۸۸٣۴ 
DT= SAMPLING INTERVAL 

NFINAL=1 

READ 1((F CIRSICO)S ICEISNOS IRS1,N) 
READ 1 (D(IO»Iz21s5N) | 
FORMAT ((8F10.0)) 

ТМ=0.0 

РЕТМТ 399 „ОТ CCFCIRS ICO) У 1С=1»М ) »з 16 21 N) 
FORMAT(///4H DTz, 1F5.3⁄////s8H Ғ(1І»49)ғ/эа(4Е16.7/)) 
۲٢۴٠٣٢٣ 399۲) E 
КРОКМАТ(///6Н О(Т)=/, 4 (1Е16.7/)} ) 

DO 400 IR-21,N 

DO 400 ICz1,N 

ТЕКМ (1В›:С) =0.0 

МОЕМ(ІК»ІС)-0.0 

ТЕКМ(ІРБ»ІК)г1.0 
TELM(IRsIC)=TERM(IRsIC)%*DT 
ОЕГР(ТЕ ,ТС ! <ТЕ МС ТК ,ТС ) 
ОЕМ (ТК ,ТС) =0.0 

ОЕ- (ТК) <0 4.0 

PHI( IRs IC}=TERM( IR IC) 

ТМ= 1. 0+ТМ 

ро 500 IR=1sN 

DO ο ο 10519 

DO 06 JN=1 >N 
DELM(IRsIC)=DELM(C IRs ICI“$TELM( IRs IN) ¥F( INS IC)¥DI/( TM+160) 
WORM(TRe IC) =TERM( IRs IN) *F (INs ITC) *¥OT/TM+WORM( IRs IC) 
DO 401 1К=1›М 

DO 401 ICz21,N 

TERMCLIRsS IC) zWORMCIRS IC) 
TELMCIRSIC)ZDELM(IR$9IC) 
DELPCIRSIC)SDELPCIRSTIC)XTELMUDLRSTC) 

РНТ (ТЕ „ТС ) <РНТ (ТВ, ТС) + ТЕКМ (ТК, ТС) 

ОЕМ ( ТК ,ТС) #0 60 

МОКМ ( ТВ ,ТС) + 0.0 

АВС-0.0 

00 21-1»М 

ОО 22/-1»М 

ΑΑΞΤΕΝΜΕΙ 9η) 

AB=ABSF (AA) 

IF(ABC-AB)3 3562 

ABC=AB 

CONTINUE 

ІЕ(0.000000005--АВС)5»5»96 

ба та 4 

РКІМТ 502» ((РНІ(ІБ»ІС)»ІСг-1» М), ІК-1.М) 
FORMAT(//9X,8HPHI(IsJ)///0((4E1627))) 
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320 
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PROGRAM MISSILE CONTINUED 


BO 800 ἐΞ}Ν 

DO 600 J=1 >N 
DFL(I)=DEL(I)+PHI(I,J)%*DELP (J yK)%D(K) 
PRINT 503 (DELCI)sT=19N) 
FORMAT(///9X» SE cA. 9)//) 
RETURN 

END 

END 


PROGRAM MISSILE CONTINUED 


EVASION 


YTDOT -ҮТООТ410.. 
EVASION LEFT-RIGHT-LEFT 
IF (KK-30) 32053215321 
YTDOT=YTDOT+20. 
GOTO 323 
ΥΤΟΟΤΞΥΤΟΟΤ-20. , 
CONTINUE 
17 (УТ2-5000.)} 32453245325 
YTDOT=YTDOT+30, | 
ЈЕ (YTZ-20000) 32693269325 
YTDOT=0.0 
CONTINUE 
END OF EVASION 
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APPENDIX II 


The missile guidance simulation is a transfer function relating 
(y) and (0). Since the objective of this study is to examine the 
effect of a digital filter, an arbitrary transfer function was se- 
lected which could represent any present day missile. This 
transfer function contains a time delay due to the radar, a fil- 
tering action due to the control (hydraulics etc.), and a second 


order system representing the missile dynamics. 


o/y = A . B i C (А-1) 
n T c DN -π--ε----- 
oS сите S“+ DS + Е 


As mentioned above the intention here is to demonstrate the 
capability of the digital filter in the presence of a large mag- 
nitude of noise. The values of the gains and time constants were 
selected for equation (A-1) after examining several different 
missile systems which use proportional navigation. Refer to 

[13 ] for further elaboration on determining the F and D matrices 


for a particular system. 
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